Hands-on Exercise 3

Description

Learning to perform Order Spatial Point Pattern Analysis in R.

Changelog

31 Aug 24: Completed Hands-on Exercise 3.

Importing Data

pacman::p_load(sf, spatstat, raster, tidyverse, tmap)

also installing the dependencies 'deldir', 'polyclip', 'spatstat.sparse', 'goftest', 'tensor', 'spatstat.data', 'spatstat.univar', 'spatstat.geom', 'spatstat.random', 'spatstat.explore', 'spatstat.model', 'spatstat.linnet', 'spatstat.utils'

Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.4:
  cannot open URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.4/PACKAGES'

package 'deldir' successfully unpacked and MD5 sums checked
package 'polyclip' successfully unpacked and MD5 sums checked
package 'spatstat.sparse' successfully unpacked and MD5 sums checked
package 'goftest' successfully unpacked and MD5 sums checked
package 'tensor' successfully unpacked and MD5 sums checked
package 'spatstat.data' successfully unpacked and MD5 sums checked
package 'spatstat.univar' successfully unpacked and MD5 sums checked
package 'spatstat.geom' successfully unpacked and MD5 sums checked
package 'spatstat.random' successfully unpacked and MD5 sums checked
package 'spatstat.explore' successfully unpacked and MD5 sums checked
package 'spatstat.model' successfully unpacked and MD5 sums checked
package 'spatstat.linnet' successfully unpacked and MD5 sums checked
package 'spatstat.utils' successfully unpacked and MD5 sums checked
package 'spatstat' successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\dexto\AppData\Local\Temp\RtmpELdYLR\downloaded_packages

spatstat installed

childcare_sf = st_read("data/child-care-services-geojson.geojson") %>%
  st_transform(crs = 3414)

Reading layer `child-care-services-geojson' from data source 
  `C:\0SMU\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

mpsz_sf = st_read(dsn = "data", layer = "MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\0SMU\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

sg_sf = st_read(dsn = "data", layer = "CostalOutline")

Reading layer `CostalOutline' from data source 
  `C:\0SMU\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data' using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21

st_crs(mpsz_sf)

Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

st_crs(sg_sf)

Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

st_crs(childcare_sf)

Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

Both mpsz_sf and sg_sf have the wrong EPSG, and must be set to 3414:

mpsz_sf = st_set_crs(mpsz_sf, 3414)

Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that

st_crs(mpsz_sf)

Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

sg_sf = st_set_crs(sg_sf, 3414)

Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that

st_crs(sg_sf)

Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

Mapping geospatial data sets

tmap_mode('view')

tmap mode set to interactive viewing

tm_shape(childcare_sf) + tm_dots()

Geospatial Data Wrangling

We need to convert the data from sf data frames into sp Spatial class:

childcare = as_Spatial(childcare_sf)
mpsz = as_Spatial(mpsz_sf)
sg = as_Spatial(sg_sf)

childcare

class       : SpatialPointsDataFrame 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :    Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           Description 
min values  :   kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>018989</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>1, MARINA BOULEVARD, #B1 - 01, ONE MARINA BOULEVARD, SINGAPORE 018989</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>THE LITTLE SKOOL-HOUSE INTERNATIONAL PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>08F73931F4A691F4</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
max values  : kml_999,                  <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>829646</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>200, PONGGOL SEVENTEENTH AVENUE, SINGAPORE 829646</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>RAFFLES KIDZ @ PUNGGOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>379D017BF244B0FA</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 

spatstat requires data in ppp object form, but there is no way to convert Spatial classes into ppp directly. As such, we must convert them into Spatial objects first:

childcare_sp = as(childcare, 'SpatialPoints')
sg_sp = as(sg, 'SpatialPolygons')

childcare_sp

class       : SpatialPoints 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

sg_sp

class       : SpatialPolygons 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

#childcare_ppp = as.ppp(childcare_sp)

As the current code gives an error, I have elected to proceed to directly take the sf dataframe and convert it into ppp directly, taking reference from https://www.paulamoraga.com/book-spatial/the-spatstat-package.html:

childcare_ppp = as.ppp(st_coordinates(childcare_sf), st_bbox(childcare_sf)) 

Warning: data contain duplicated points

childcare_ppp

Marked planar point pattern: 1545 points
marks are numeric, of storage type  'double'
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units

plot(childcare_ppp)

summary(childcare_ppp)

Marked planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 11 decimal places

marks are numeric, of type 'double'
Summary:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0 

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units

Handling duplicated points

Check for duplication in ppp object:

any(duplicated(childcare_ppp))

[1] TRUE

We can check for the number of points using multiplicity():

sum(multiplicity(childcare_ppp) > 1)

[1] 128

Plotting this out, we can identify the duplicates by the colour of the dots. The darker dots signals multiple points with the same coordinates:

tm_shape(childcare) + tm_dots(alpha=0.3)

There are 3 techniques to deal with duplicates:

Delete duplicates - May lose some important data

Jittering - Add small perturbation to duplicate points to make them not occupy the exact same space

Marks - Make each point unique and attach duplicates as marks (attributes of the points)

I have applied the jittering approach:

childcare_ppp_jit = rjitter(childcare_ppp, retry = TRUE, nsim = 1, drop = TRUE)
any(duplicated(childcare_ppp_jit))

[1] FALSE

Next, we will convert sg SpatialPolygon object into owin:

sg_owin = as.owin(sg_sf)
plot(sg_owin)

Now, we combine the owin map with the childcare data:

childcareSG_ppp = childcare_ppp[sg_owin]
summary(childcareSG_ppp)

Marked planar point pattern:  1545 points
Average intensity 2.129929e-06 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 11 decimal places

marks are numeric, of type 'double'
Summary:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0 

Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401

Kernel Density Estimation (KDE)

Usage of automatic bandwidth selection method to compute kernel density:

kde_childcareSG_bw = density(childcareSG_ppp, sigma = bw.diggle, edge = TRUE,
                             kernel = "gaussian")
plot(kde_childcareSG_bw)

Retrieving the bandwidth used:

bw = bw.diggle(childcareSG_ppp)
bw

   sigma 
298.4095 

Rescalling KDE values to change meters to kilometers:

childcareSG_ppp.km = rescale.ppp(childcareSG_ppp, 1000, "km")
kde_childcareSG_bw.km = density(childcareSG_ppp.km, sigma = bw.diggle, edge = TRUE,
                                kernel = "gaussian")
plot(kde_childcareSG_bw.km)

Using different bandwidth calculation methods:

bw.CvL(childcareSG_ppp.km)

   sigma 
4.543278 

bw.scott(childcareSG_ppp.km)

 sigma.x  sigma.y 
2.224898 1.450966 

bw.ppl(childcareSG_ppp.km)

    sigma 
0.3897114 

bw.diggle(childcareSG_ppp.km)

    sigma 
0.2984095 

Comparing diggle vs ppl

kde_childcareSG_ppl.km = density(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE,
                                kernel = "gaussian")
par(mfrow=c(1,2), mai=c(0.2,0.2,0.2,0.2))
plot(kde_childcareSG_bw.km, main = "diggle")
plot(kde_childcareSG_ppl.km, main = "ppl")

Using different kernel methods:

par(mfrow=c(2,2), mai=c(0.2,0.2,0.2,0.2))

plot(density(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE,
                                kernel = "gaussian"), main = "Gaussian")
plot(density(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE,
                                kernel = "epanechnikov"), main = "Epanechnikov")

Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel

plot(density(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE,
                                kernel = "quartic"), main = "Quartic")

Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel

plot(density(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE,
                                kernel = "disc"), main = "Disc")

Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel

Computing KDE using fixed bandwidth:

kde_childcareSG_600 = density(childcareSG_ppp.km, sigma = 0.6, edge = TRUE, 
                                 kernel = "gaussian")
plot(kde_childcareSG_600)

Computing KDE using adaptive bandwidth:

kde_childcareSG_adaptive = adaptive.density(childcareSG_ppp.km, method = "kernel")
plot(kde_childcareSG_adaptive)

Now we convert the object into a grid object for mapping:

gridded_kde_childcareSG_bw <- as(kde_childcareSG_bw.km, "SpatialGridDataFrame")
spplot(gridded_kde_childcareSG_bw)

Convert gridded kernal density objects into RasterLayer object:

kde_childcareSG_bw_raster = raster(gridded_kde_childcareSG_bw)
kde_childcareSG_bw_raster

class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -8.476185e-15, 28.51831  (min, max)

Assign a projection system to RasterLayer, as our current object has no crs:

projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster

class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : v 
values     : -8.476185e-15, 28.51831  (min, max)

tmap_mode("plot")

tmap mode set to plotting

tm_shape(kde_childcareSG_bw_raster) + 
  tm_raster("v", palette = "viridis") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

pg = mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm = mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck = mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw = mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")
plot(pg, main = "Ponggol")

Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(tm, main = "Tampines")

Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(ck, main = "Choa Chu Kang")

Warning: plotting the first 10 out of 15 attributes; use max.plot = 15 to plot
all

plot(jw, main = "Jurong West")

Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

Converting the sf objects into owin objects and combining it with childcare points, including the rescalling to km:

pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)
childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]
childcare_pg_ppp.km = rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale.ppp(childcare_jw_ppp, 1000, "km")
par(mfrow=c(2,2),mai=c(0.2,0.2,0.2,0.2))
plot(childcare_pg_ppp.km, main="Punggol")
plot(childcare_tm_ppp.km, main="Tampines")
plot(childcare_ck_ppp.km, main="Choa Chu Kang")
plot(childcare_jw_ppp.km, main="Jurong West")

Plotting the KDE with diggle for the 4 areas:

par(mfrow=c(2,2),mai=c(0.2,0.2,0.2,0.2))
plot(density(childcare_pg_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tempines")
plot(density(childcare_ck_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Jurong West")

Plotting the KDE with fixed bandwidth for the 4 areas:

par(mfrow=c(2,2),mai=c(0.2,0.2,0.2,0.2))
plot(density(childcare_pg_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Chou Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Jurong West")
plot(density(childcare_pg_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tampines")

Nearest Neighbour Analysis

We are using the Clark-Evans test <clarkevans.test()> to perform a test on the data. The hypotheses are as follows:

Ho = The distribution of childcare services are randomly distributed.

H1= The distribution of childcare services are not randomly distributed.

The 95% confident interval will be used.

Testing spatial point patterns

clarkevans.test(childcareSG_ppp, correction = "none", clipregion = "sg_owin", 
                alternative = c("clustered"), nsim = 99)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.55631, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

CCK planning area

clarkevans.test(childcare_ck_ppp, correction = "none", clipregion = NULL, 
                alternative = c("two.sided"), nsim = 999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.93706, p-value = 0.347
alternative hypothesis: two-sided

Tampines planning area

clarkevans.test(childcare_tm_ppp, correction = "none", clipregion = NULL, 
                alternative = c("two.sided"), nsim = 999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.80594, p-value = 0.0004612
alternative hypothesis: two-sided

G-Function

Computing G-function using Gest():

#CCK
G_CK = Gest(childcare_ck_ppp, correction = "border")
plot(G_CK, xlim = c(0,500))

#Tampines
G_TM = Gest(childcare_tm_ppp, correction = "best")
plot(G_TM, xlim = c(0,500))

Spatial Randomness Test

Usage of the Monte Carlo test to test the random distribution of childcare services in CCK and Tampines.

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

G_CK.csr = envelope(childcare_ck_ppp, Gest, nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(G_CK.csr, xlim = c(0,500))

G_TM.csr = envelope(childcare_tm_ppp, Gest, correction = "all", nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(G_TM.csr, xlim = c(0,500))

F-Function

Compute F-function using Fest():

F_CK = Fest(childcare_ck_ppp)
plot(F_CK)

F_TM = Fest(childcare_tm_ppp, correction = "best")
plot(F_TM)

Spatial Randomness Test

Usage of the same hypothesis as shown in G-Function’s test.

F_CK.csr = envelope(childcare_ck_ppp, Fest, nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(F_CK.csr)

F_TM.csr = envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(F_TM.csr)

K-Function

Compute K-function using Kest():

K_CK = Kest(childcare_ck_ppp, correction = "Ripley") 
plot(K_CK, . - r ~ r, ylab = "K(d)-r", xlab = "d(m)")

K_TM = Kest(childcare_tm_ppp, correction = "Ripley") 
plot(K_TM, . - r ~ r, ylab = "K(d)-r", xlab = "d(m)", xlim = c(0,1000))

Spatial Randomness Test

Usage of the same hypothesis as shown in G-Function’s test.

K_CK.csr = envelope(childcare_ck_ppp, Kest, nsim = 99, rank = 1, glocal = TRUE) 

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(K_CK.csr, . - r ~ r, xlab = "d", ylab = "K(d)-r")

K_TM.csr = envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal = TRUE) 

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(K_TM.csr, . - r ~ r, xlab = "d", ylab = "K(d)-r", xlim = c(0,500))

L-Function

Compute L-function using Lest():

L_CK = Lest(childcare_ck_ppp, correction = "Ripley") 
plot(L_CK, . - r ~ r, ylab = "L(d)-r", xlab = "d(m)")

L_TM = Lest(childcare_tm_ppp, correction = "Ripley") 
plot(L_TM, . - r ~ r, ylab = "L(d)-r", xlab = "d(m)", xlim = c(0,1000))

Spatial Randomness Test

Usage of the same hypothesis as shown in G-Function’s test.

L_CK.csr = envelope(childcare_ck_ppp, Lest, nsim = 99, rank = 1, glocal = TRUE) 

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(L_CK.csr, . - r ~ r, xlab = "d", ylab = "L(d)-r")

L_TM.csr = envelope(childcare_tm_ppp, Lest, nsim = 99, rank = 1, glocal = TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(L_TM.csr, . - r ~ r, xlab = "d", ylab = "L(d)-r", xlim = c(0,500))